#!git clone https://github.com/TPM034A/Q2_2025
#!pip install -r Q2_2025/requirements_colab.txt
#!mv "/content/Q2_2025/Lab_sessions/lab_session_03/data" /content/data

## IMPORTANT: You have to be on the TUDelft network (eduroam) or under eduVPN to run this script
download_full_dataset = False

from assets import image_downloader as imd
imd.download_images(download_full_dataset)

Downloading images...
Download complete!
Unzipping images...
Unzip complete!
Removing zip file...
Done!

# Basic libraries
import numpy as np
import pandas as pd
import geopandas as gpd

# ML tools
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from sklearn.model_selection import train_test_split, GridSearchCV, cross_validate
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error, make_scorer,log_loss

# Models
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor, VotingRegressor
from sklearn.neural_network import MLPRegressor

# Visualization libraries
import seaborn as sns
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import ImageGrid
from matplotlib.ticker import FixedLocator
from sklearn.tree import plot_tree
from branca.element import Figure

# Other libraries
from pathlib import Path
from shapely.geometry import Point
from PIL import Image
import pickle

# Data folder path
data = Path('data')

# Reading image df
img_metadata = pd.read_csv(data/'image_tabular'/'image_metadata.csv')

# Reading hexagons gdf
hexagons = gpd.read_file(data/'geo'/'hexagons.gpkg')

# Verify geographic coordinate system / projection
print(hexagons.crs)

EPSG:4326

img_metadata.head(5)

img_metadata.describe()

# Histogram of the attractiveness at image level 
fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(img_metadata['attractiveness'], bins=100)
ax.set_xlabel('Visual attractiveness')
ax.set_ylabel('Frequency')
ax.set_title('Histogram of the attractiveness at image level')
plt.show()

def dataframe_to_geodataframe_nl(original_dataframe, latitud_column_name, longitude_column_name):
    '''
    This Function converts a dataframe into a geodataframe using the latitud and longitude columns.
    The output will be ready to use the plot function from geopandas.
    '''

    ## Creating the Point geometry using lat/lng columns
    original_dataframe['geometry'] = original_dataframe.apply(lambda x: Point(x[longitude_column_name], x[latitud_column_name]), axis=1)

    ## Creating the geodataframe (we used crs 4326 because it is the code for the latitud and longitude)
    geodataframe = gpd.GeoDataFrame(original_dataframe, geometry='geometry', crs=4326)

    ## Changing the crs to the same as the hexagons This is 28992, the projection used in the Netherlands.
    geodataframe = geodataframe.to_crs(28992)

    return geodataframe

# ADD HERE YOUR ANSWER TO EXERCISE 1

# A. Checking the number of images per hexagon
hex_num_imgs = img_metadata[['hex_id', 'img_id']].groupby('hex_id').count()
hex_num_imgs = hex_num_imgs.rename(columns={'img_id': 'num_imgs'})

fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(hex_num_imgs['num_imgs'])
ax.set_xlabel('Number of images per hexagon')
ax.set_ylabel('Frequency')
ax.set_title('Histogram of the number of images per hexagon')
plt.show()

# B. Creating the histogram of images by year and month
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
ax[0].hist(img_metadata['year'], bins=15)
ax[0].set_xlabel('Year')
ax[0].set_ylabel('Frequency')
ax[0].set_title('Histogram of the year of the images')
ax[1].hist(img_metadata['month'], bins=12)
ax[1].set_xlabel('Month')
ax[1].set_ylabel('Frequency')
ax[1].set_title('Histogram of the month of the images')
plt.show()

# The histogram by years shows a gap between 2011 and 2013. Also 2014 is the year with more images. 
# The histogram by month shows that most images are from summer, especially from July.

# C. To visualize the spatial distribution of the data, we can convert the dataframe into a geodataframe using the function provided in the lab session.

img_metadata_gdf = dataframe_to_geodataframe_nl(img_metadata, 'lat', 'lng')

## Now we can plot the data

# For years
fig, axes = plt.subplots(nrows=3, ncols=5, figsize=(18, 10))
for ax, year in zip(axes.flatten(), range(2008, 2023)):
    if len(img_metadata[img_metadata.year == year]) > 0:
        img_metadata_gdf[img_metadata.year == year].plot(ax=ax, markersize=0.5, alpha=0.5)
    ax.set_title('Year {}'.format(year))
    ax.set_axis_off()
plt.show()

# The same for each month
fig, axes = plt.subplots(nrows=3, ncols=4, figsize=(15, 10))
for ax, month in zip(axes.flatten(), range(1, 13)):
    if len([img_metadata.month == month]) > 0:
        img_metadata_gdf[img_metadata.month == month].plot(ax=ax, markersize=0.5, alpha=0.5)
    ax.set_title('Month: {}'.format(month))
    ax.set_axis_off()
plt.show()

# The maps show that from 2014 images have been taken broadly across the city. Most of the images are from summer time where also there is a broader coverage.

# D. Summer pictures maybe are more attractive than winter pictures beacuse there are more green vegetation.

hexagons.head(5)

fig = Figure(width=400, height=300)
fig.add_child(hexagons.explore(marker_kwds={'radius': 10}, zoom_start=15))
display(fig)

# ADD HERE YOUR ANSWER TO EXERCISE 1

# A. Getting the average attractiveness of each hexagon
hex_avg_attrac = img_metadata[['hex_id', 'attractiveness']].groupby('hex_id').mean()

# B. Adding the average attractiveness levels to the hexagons of the  geodataframe
hex_avg_attrac = pd.merge(hexagons, hex_avg_attrac, on='hex_id', how='left')

# Now we have an average attractiveness for each hexagon.
hex_avg_attrac.head(5)

# C. Plotting the average attractiveness of hexagons
fig = Figure(width=600, height=400)

fig.add_child(hex_avg_attrac.explore(column = 'attractiveness', cmap = 'RdYlGn', scheme = 'quantiles'))
display(fig)

# D. There is a clear pattern of attractiveness in the city. The most relatively unattractive areas are the city center compared to the periphery. 
# The most attractive areas are the city's periphery, especially the south and the northeast, where parks and green areas are located.

image_path = lambda file: data/'images'/file
visual_images = img_metadata.loc[img_metadata['in_folder'] == 1].copy()

# Create a figure
fig, ax = plt.subplots(1,2,figsize=(12,3.5))

# Draw a random image
rnd_img = visual_images.sample(1)[['img_path','attractiveness']]

# Visualising the random image
img = Image.open(image_path(rnd_img.iloc[0]['img_path'])).convert('RGB')
ax[0].imshow(img)

# Remove the axis
# ax[0].axis('off')
ax[0].set_xlabel('Random image')
ax[0].set_xticks([])
ax[0].set_yticks([])

# Add the histogram of the attractiveness of all 101k images
sns.histplot(img_metadata['attractiveness'],bins=50, stat='density', ax=ax[1])
ax[1].set_xlabel('Visual attractiveness')

# Show where its attractiveness is located in the distribution
img_attractiveness = rnd_img.iloc[0]['attractiveness']
bin_edges = ax[1].patches
for patch in bin_edges:
    if patch.get_x() <= img_attractiveness <= patch.get_x() + patch.get_width():
        patch.set_facecolor('red')
plt.show()

# Set the number of percentiles and images per row
n_percentiles = 6
images_per_row = 5

# Create a new column with the quantile
visual_images['quantile'], quantile_boundaries = pd.qcut(visual_images['attractiveness'], n_percentiles, labels=False, retbins=True)

# Create a figure with a grid
fig = plt.figure(figsize=(20., 100.))
grid = ImageGrid(fig, 111, nrows_ncols=(n_percentiles,images_per_row+1),axes_pad=0.05)

# 
for i in range(n_percentiles):
    # Draw 'images_per_row' random images from the quantile i
    auxs = visual_images[visual_images['quantile']==i].sample(images_per_row)[['img_path', 'attractiveness']].reset_index(drop=True)
    
    # Plot the images in the grid
    for j, row in auxs.iterrows():
        img = Image.open(image_path(row['img_path'])).convert('RGB')
        grid[(i*(images_per_row+1))+j+1].imshow(img)
        grid[(i*(images_per_row+1))+j+1].axis('off')

    # Add the quantile label and the quantile boundaries     
    grid[i*(images_per_row+1)].text(-300, 300, f'Quantile {i+1}\n ({quantile_boundaries[i]:0.2f}, {quantile_boundaries[i+1]:0.2f})', fontsize=15, ha='center', va='center')
    
    # Do not show axis
    grid[i*(images_per_row+1)].axis('off')

# Show the figure
plt.show()

# Set the number of images per row
images_per_row = 5

# Sort the images by attractiveness and get the first and last images
mins = visual_images.sort_values(by='attractiveness', ascending=True).head(images_per_row*2).reset_index(drop=True)

# Create a figure with a grid
fig = plt.figure(figsize=(20., 100.))
grid = ImageGrid(fig, 111, nrows_ncols=(2,images_per_row),axes_pad=0.05)

# Plot the images in the grid
for i, row in mins.iterrows():
    img = Image.open(image_path(row['img_path'])).convert('RGB')
    grid[i].imshow(img)
    
    # Do not show axis
    grid[i].axis('off')

# Show the figure
plt.show()

# Set the number of images per row
images_per_row = 5

# Sort the images by attractiveness and get the first and last images
maxs = visual_images.sort_values(by='attractiveness', ascending=True).tail(images_per_row*2).reset_index(drop=True)

# Create a figure with a grid
fig = plt.figure(figsize=(20., 100.))
grid = ImageGrid(fig, 111, nrows_ncols=(2,images_per_row),axes_pad=0.05)

# Plot the images in the grid
for i, row in maxs.iterrows():
    img = Image.open(image_path(row['img_path'])).convert('RGB')
    grid[i].imshow(img)
    
    # Do not show axis
    grid[i].axis('off')

# Show the figure
plt.show()

# Load the embeddings
embeddings = pd.read_csv(data/'image_tabular'/'image_embeddings.csv')

embeddings.head()

embeddings.describe()

embeddings_ym = pd.merge(embeddings, img_metadata[['img_id','year','month', 'attractiveness']], on='img_id')

# heatmap of correlations
# Create plot
fig, axes = plt.subplots(figsize=(7, 6))
fig.set_tight_layout(True)

# Compute correlation matrix
corr = embeddings_ym.drop(columns = 'img_id').corr()

# Create upper triangular matrix to mask the upper triangular part of the heatmap
corr_mask = np.triu(np.ones_like(corr, dtype=bool))

# Generate a custom diverging colormap (because it looks better)
corr_cmap = sns.diverging_palette(230, 20, as_cmap=True)

sns.heatmap(corr, mask = corr_mask, cmap=corr_cmap, annot=False,square = True, linewidths=.5, ax = axes)
plt.show()

# In the case of attractiveness we can see that the correlation is in some way correlated with some features. The case for year and month is almost close to zero.

# Merging the embeddings with the original dataset
image_embeddings = pd.merge(img_metadata, embeddings, on='img_id')
image_embeddings.head(5)

# Generating the embedding columns
df = image_embeddings.copy()
emb_cols = embeddings.drop(columns=['img_id']).columns
features = emb_cols.to_list()
predict = 'attractiveness'

# Setting all dataset between 0 and 1
scaler = StandardScaler()
df[features + [predict]] = scaler.fit_transform(df[features + [predict]])

# Splitting the dataset
X_train, X_test, Y_train, Y_test = train_test_split(df[features], df[predict], test_size=0.2, random_state=44)

def eval_regression_perf(model, X_train, X_test, Y_train, Y_test):
    
    # Make prediction with the trained model
    Y_pred_train = model.predict(X_train)
    Y_pred_test = model.predict(X_test)

    # Create a function that computes the MSE, MAE, and R2
    def perfs(Y,Y_pred):
        mse = mean_squared_error(Y,Y_pred)
        mae = mean_absolute_error(Y,Y_pred)
        R2 = r2_score(Y,Y_pred)
        return mse,mae,R2

    # Apply the perfs function to the train and test data sets
    mse_train, mae_train, r2_train = perfs(Y_train,Y_pred_train)
    mse_test,  mae_test , r2_test  = perfs(Y_test,Y_pred_test)
        
    # Print results
    print('Performance')
    print(f'Mean Squared  Error Train | Test: \t{mse_train:>7.4f}\t|  {mse_test:>7.4f}')
    print(f'Mean Absolute Error Train | Test: \t{mae_train:>7.4f}\t|  {mae_test:>7.4f}')
    print(f'R2                  Train | Test: \t{ r2_train:>7.4f}\t|  {r2_test:>7.4f}\n')

# Create and train a multiple linear regression model
regr = LinearRegression(fit_intercept = True)

# Fit the model on the training data
regr.fit(X_train,Y_train)

# Evaluate the performance of the trained model
eval_regression_perf(regr, X_train, X_test, Y_train, Y_test)

# We take an ML approach (i.e. we use a train and tests set). But, we can still look at the regression coefficients to see which features are more important than others.
regr_results_int = pd.DataFrame(index = ['Intercept'], data = regr.intercept_, columns = ['Coefficient'])
regr_results_wei = pd.DataFrame(index = df[features].columns, data = regr.coef_, columns = ['Coefficient'])
regr_results = pd.concat([regr_results_int,regr_results_wei],axis=0)

# Show the results
regr_results.sort_values(by='Coefficient', ascending=False).plot(kind='bar', figsize=(10, 4))

Performance
Mean Squared  Error Train | Test: 	 0.3485	|   0.3515
Mean Absolute Error Train | Test: 	 0.4624	|   0.4621
R2                  Train | Test: 	 0.6504	|   0.6528

<Axes: >

# Create RF object for hyperparameter tuning
rf_gs = RandomForestRegressor(n_estimators=50,criterion= "squared_error",random_state=5)

# Define the hyperparameter search space
hyperparameter_space = {
    'max_depth': [5, 20], 
    'max_leaf_nodes': [500,2500],
    # 'min_samples_leaf': [15],     # Not tuned here too avoid long computation time 
    # 'max_features':[0.7]          # Not tuned here too avoid long computation time
    } 

# Create scoring function
scorer = make_scorer(mean_squared_error, greater_is_better = False)

# Create the grid_search object, with using the MLP classifier
folds = 5 # Number of cross validation splits
rf_gridsearch = GridSearchCV(rf_gs, hyperparameter_space, n_jobs=-1, cv=folds, scoring=scorer, return_train_score=True)

# if the model RF_gridsearch.pickle does not exist, conduct the gridsearch, else load the grid_search object
if not Path('RF_gridsearch.pickle').exists():
    
    # Fit the grid search object to the training data
    # Note that this is computationally expensive! 
    # It may take up to 5 minutes, since multiple models need to be trained multiple times
    rf_gridsearch.fit(X_train, Y_train)

    # Save the trained grid_search object
    with open('RF_gridsearch.pickle', 'wb') as f:
        pickle.dump(rf_gridsearch, f)
else:
    # Load the trained grid_search object
    with open('RF_gridsearch.pickle', 'rb') as f:
        rf_gridsearch = pickle.load(f)

# Print the best parameters
print('Best parameters:')
print(f'Tuned hyperparameters:\n{rf_gridsearch.best_params_}')
print(f'Mean Squared Error:\t{(-rf_gridsearch.best_score_):0.3f}')

Best parameters:
Tuned hyperparameters:
{'max_depth': 20, 'max_leaf_nodes': 2500}
Mean Squared Error:	0.373

# Create a dataframe from the results
df_results = pd.DataFrame(rf_gridsearch.cv_results_)

# Extract relevant columns
df_results = df_results[['params', 'mean_test_score', 'mean_train_score']]

# Display the top performing models in a plot
top_models = df_results.sort_values(by='mean_test_score', ascending=False)

# Create a figure
fig, ax = plt.subplots(figsize=(10, 3))

# Convert the dict row names to a list of strings, and replace the commas with newlines
param_list = [str(dict_) for dict_ in top_models['params']]
param_list = [param.replace(',', '\n') for param in param_list]
param_list = [param.replace('{', '')   for param in param_list]
param_list = [param.replace('}', '')   for param in param_list]

# Create a barplot
sns.barplot(x= param_list, y=top_models['mean_test_score'], ax=ax)

# Set the x-axis labels to the parameter list
ax.set_xlabel('Hyperparameters')

ax.xaxis.set_major_locator(FixedLocator(range(len(param_list))))
ax.set_xticklabels(param_list, rotation=90)

# Show the plot
plt.show()

# Set the hyperparameters of the Random Forest to the best values found, e.g. 
max_depth = rf_gridsearch.best_params_['max_depth']
max_leaf_nodes = rf_gridsearch.best_params_['max_leaf_nodes']
max_features = 0.7
min_samples_leaf = 15

# Create the Random Forest object with the best hyperparameters
rf = RandomForestRegressor(n_estimators=100, max_depth=max_depth, max_features=max_features, max_leaf_nodes=max_leaf_nodes, min_samples_leaf=min_samples_leaf, random_state=0,n_jobs=-1)

# Train the Randon Forest
rf.fit(X_train,Y_train)

# Evaluate the performance of the hyperparameter optimised RF model
eval_regression_perf(rf,X_train,X_test, Y_train, Y_test)

Performance
Mean Squared  Error Train | Test: 	 0.2385	|   0.3714
Mean Absolute Error Train | Test: 	 0.3789	|   0.4762
R2                  Train | Test: 	 0.7607	|   0.6332

# Plot the feature importances
fig, ax = plt.subplots(figsize=(8,10))
sns.barplot(y=features,x=rf.feature_importances_, ax = ax).set(title='Feature importance to predict Attractiveness')
ax.set_xlabel('Feature importance')

# Print most important features
sorted_indices = np.argsort(rf.feature_importances_)[::-1]
most_imp_features_rf = [features[i] for i in sorted_indices]
print(f'Top 5 most important features:\n {most_imp_features_rf[:6]}')

# Show plot
plt.show()

Top 5 most important features:
 ['Fhat004', 'Fhat002', 'Fhat006', 'Fhat001', 'Fhat009', 'Fhat014']

# Create MLP object (plain vanilla MLP)
mlp_gs = MLPRegressor(activation = 'tanh', solver='adam', batch_size=250, max_iter=2000)

# Define the hyperparameter search space
# 'hidden_layer_sizes' defines the number of nodes and layers
# 'alpha' governs the L2 regularisation
# 'learning_rate_init' governs the learning rate.
hyperparameter_space = {
    
    'hidden_layer_sizes': [(20,20),(30,30),(25,25)],
    'alpha': [1,0.1,0.01],
    'learning_rate_init': [0.01, 0.001,0.0001, 0.00001]}

# Create scoring function
scorer = make_scorer(mean_squared_error, greater_is_better = False)

# Create the grid_search object, with using the MLP classifier
folds = 5 # Number of cross validation splits
mlp_gridsearch = GridSearchCV(mlp_gs, hyperparameter_space, n_jobs=-1, cv=folds,scoring = scorer)

# if the model mlp_gridsearch.pickle does not exist, conduct the gridsearch, else load the grid_search object
if not Path('mlp_gridsearch.pickle').exists():
    
    # Fit the grid search object to the training data
    # Note that this is computationally expensive! 
    # It may take up to 5 minutes, since multiple models need to be trained multiple times
    mlp_gridsearch.fit(X_train, Y_train)

    # Save the trained grid_search object
    with open('mlp_gridsearch.pickle', 'wb') as f:
        pickle.dump(mlp_gridsearch, f)
else:
    # Load the trained grid_search object
    with open('mlp_gridsearch.pickle', 'rb') as f:
        mlp_gridsearch = pickle.load(f)

print(f'Optimal hyperparameters:\n{mlp_gridsearch.best_params_}')
print(f'Mean Squared Error:\t{(-mlp_gridsearch.best_score_):0.3f}')

Optimal hyperparameters:
{'alpha': 0.1, 'hidden_layer_sizes': (30, 30), 'learning_rate_init': 0.001}
Mean Squared Error:	0.289

# Create and train a MLP model (sklearn)
hidden_layer_sizes = mlp_gridsearch.best_params_['hidden_layer_sizes']
alpha = mlp_gridsearch.best_params_['alpha']
learning_rate_init = mlp_gridsearch.best_params_['learning_rate_init']

mlp_gs = MLPRegressor(hidden_layer_sizes = hidden_layer_sizes, alpha = alpha, learning_rate_init = learning_rate_init, activation='tanh', solver='adam', batch_size = 250 , max_iter=2000, random_state=1)
mlp_gs.fit(X_train,Y_train.to_numpy().ravel())

# Evaluate the performance of the trained model
eval_regression_perf(mlp_gs, X_train, X_test, Y_train.to_numpy().ravel(), Y_test.to_numpy().ravel())

Performance
Mean Squared  Error Train | Test: 	 0.2801	|   0.2896
Mean Absolute Error Train | Test: 	 0.4127	|   0.4187
R2                  Train | Test: 	 0.7190	|   0.7140

# Scatter true shares (x-axis) against the predicted shares by the different ML models (y-axis)
fig, ax = plt.subplots(1,3,figsize=(18,4), sharex=True, sharey=True)
ax[0].scatter(x = Y_test, y = regr.predict(X_test),  marker = '.', s = 25, alpha = 0.2, c = '#ff7f0e', label = f'Regression | R2 = {regr.score(X_test,Y_test):0.3f}')
ax[1].scatter(x = Y_test, y = rf.predict(X_test),    marker = '.', s = 25, alpha = 0.2, c = '#2ca02c', label = f'Random Forest| R2 = {rf.score(X_test,Y_test):0.3f}')
ax[2].scatter(x = Y_test, y = mlp_gs.predict(X_test),marker = '.', s = 25, alpha = 0.2, c = '#d62728', label = f'MLP | R2 = {mlp_gs.score(X_test,Y_test):0.3f}')

# Add labels, legend, and title to each plot
for n in range(0,3,1):
    ax[n].set_xlabel('True Attractiveness')
    ax[n].set_ylabel('Predicted attractiveness')
    ax[n].set_xlim(-4,4)
    ax[n].set_ylim(-4,4)
    ax[n].legend()
    ax[n].set_title('True versus predicted attractiveness (test data)')

plt.show()

# CODE YOUR ANSWERS HERE (Use as many cells as you need)

# A: The most striking difference is that the MLP (Red) handles extreme values much better than the Regression or 
# Random Forest models. In the Random Forest and regression plots, the points flatten out at the ends (top-right and bottom-left) 
# compared to the MLP plot. This means these models fail to predict "very high" or "very low" attractiveness accurately, 
# tending instead to predict values closer to the average. The MLP plot shows a tighter diagonal cloud that extends further into
# the corners. It captures the non-linear complexity of the data better, resulting in the highest accuracy score (R2=0.7).

# B: An ideal scatter plot would show a perfect straight diagonal line rising at a 45-degree angle from the bottom-left to the 
# top-right. There would be no "cloud" or dispersion; every single point would lie exactly on that line.

# Plot kernel densities for the test data
fig, ax = plt.subplots(figsize=(10, 7))

# Plot kernel density of the true distribution, i.e of Y_test
sns.kdeplot(data = Y_test, ax=ax, legend=True,fill=True,alpha=.3, clip =(-4,4), label = 'True')

# Plot kernel densities of the predicted distributions, based on the three models
sns.kdeplot(data = regr.predict(X_test),ax=ax, legend=True,fill=True,alpha=.3, clip =(-4,4), label = f'Regression | R2 = {regr.score(X_test,Y_test):0.3f}')
sns.kdeplot(data = rf.predict(X_test),  ax=ax, legend=True,fill=True,alpha=.3, clip =(-4,4), label = f'Random Forest | R2 = {rf.score(X_test,Y_test):0.3f}')
sns.kdeplot(data = mlp_gs.predict(X_test),  ax=ax, legend=True,fill=True,alpha=.3, clip =(-4,4), label = f'MLP | R2 = {mlp_gs.score(X_test,Y_test):0.3f}')

# Add legend and title to the plot
#ax.legend(['True','Regression model', 'Decision tree','Random Forest'])
ax.legend()
plt.title('True vs predicted Attractiveness')
ax.set_xlim(-4,4)
plt.show()

# CODE YOUR ANSWERS HERE (Use as many cells as you need)
## THERE NO UNIQUE ANSWERS, BUT THIS IS AN EXAMPLE OF A GOOD ANSWER:

#A: The MLP (Multi-Layer Perceptron) is the best model for this task based on the provided scatter plots and performance metrics. 
# It achieved the highest R2 of 0.714, compared to Regression (0.653) and Random Forest (0.633). The scatter plots showed that 
# the MLP was better at capturing the non-linear relationships in the data, maintaining a diagonal shape even at the 
# extreme ends of the attractiveness scale where the other models did not.

#B: This indicates that the models are predicting the expected (average) attractiveness based on features, 
# which inherently eliminates the natural noise and variability found in individual human ratings. It tells us 
# the models are "conservative" as they fail to capture the full spectrum of extreme reactions (very ugly or 
# very beautiful places), resulting in a compressed range of predicted values compared to reality.

#C:It provided a necessary benchmark to understand the data's complexity. By comparing different models, we demonstrated 
# that the relationship between urban features and attractiveness is non-linear, which the MLP captured best. Furthermore, 
# consistent results across models increase our confidence; had we obtained contradictory trends, the validity of our findings 
# would be much weaker.

	img_id	img_path	year	month	lat	lng	hex_id	attractiveness
0	0	image_261669_s_a.png	2019	5	51.992875	4.595223	8b19694da002fff	0.745663
1	1	image_261609_s_a.png	2021	1	51.992723	4.595357	8b19694da003fff	0.729360
2	2	image_261683_s_a.png	2017	7	51.993201	4.594991	8b19694da006fff	0.767513
3	3	image_261590_s_a.png	2016	7	51.991961	4.595907	8b19694da008fff	0.721686
4	4	image_261599_s_a.png	2017	7	51.992304	4.595650	8b19694da00efff	0.572475

	img_id	year	month	lat	lng	attractiveness	in_folder
count	101444.000000	101444.000000	101444.000000	101444.000000	101444.000000	101444.000000	101444.000000
mean	50721.500000	2015.483341	6.995032	51.918793	4.490640	-0.285988	0.009858
std	29284.504691	4.358636	2.427367	0.030294	0.048692	0.529770	0.098796
min	0.000000	2008.000000	1.000000	51.852995	4.380001	-2.526708	0.000000
25%	25360.750000	2014.000000	6.000000	51.891868	4.457293	-0.669447	0.000000
50%	50721.500000	2016.000000	6.000000	51.919943	4.489459	-0.319260	0.000000
75%	76082.250000	2019.000000	9.000000	51.943226	4.530447	0.078983	0.000000
max	101443.000000	2022.000000	12.000000	51.994166	4.601524	1.748261	1.000000

	hex_id	geometry
0	8b19694da002fff	POLYGON ((4.59455 51.99293, 4.5946 51.9927, 4....
1	8b19694da003fff	POLYGON ((4.59496 51.99263, 4.595 51.9924, 4.5...
2	8b19694da006fff	POLYGON ((4.59482 51.99331, 4.59486 51.99309, ...
3	8b19694da008fff	POLYGON ((4.59577 51.99204, 4.59582 51.99181, ...
4	8b19694da00efff	POLYGON ((4.59537 51.99234, 4.59541 51.99211, ...

	hex_id	geometry	attractiveness
0	8b19694da002fff	POLYGON ((4.59455 51.99293, 4.5946 51.9927, 4....	0.270133
1	8b19694da003fff	POLYGON ((4.59496 51.99263, 4.595 51.9924, 4.5...	0.351845
2	8b19694da006fff	POLYGON ((4.59482 51.99331, 4.59486 51.99309, ...	0.791481
3	8b19694da008fff	POLYGON ((4.59577 51.99204, 4.59582 51.99181, ...	0.766068
4	8b19694da00efff	POLYGON ((4.59537 51.99234, 4.59541 51.99211, ...	0.681892

	img_id	Fhat000	Fhat001	Fhat002	Fhat003	Fhat004	Fhat005	Fhat006	Fhat007	Fhat008	...	Fhat015	Fhat016	Fhat017	Fhat018	Fhat019	Fhat020	Fhat021	Fhat022	Fhat023	Fhat024
0	0	-0.721522	0.232802	0.139053	0.332571	-0.643862	0.063346	-0.143065	0.197804	0.627572	...	0.093596	-0.024934	-0.183208	-0.329915	0.334317	-0.034518	-0.005868	-0.170634	-0.033538	-0.085577
1	1	-0.623503	0.638402	0.658155	0.192114	-0.627638	0.011822	-0.144213	0.119007	0.454119	...	0.278526	0.025043	-0.273027	-0.256410	0.200388	-0.166567	0.180788	-0.246929	0.179483	0.129283
2	2	-0.100401	0.421808	1.200465	-0.040461	-0.798148	-0.101276	-0.656213	0.018648	-0.331052	...	-0.053597	-0.061117	0.003480	-0.470245	-0.163239	0.103589	-0.206567	-0.512829	-0.079979	-0.046585
3	3	0.161000	0.808097	0.024000	0.184818	-0.909536	0.111023	-0.225173	0.245244	-0.088505	...	0.083257	-0.193682	-0.133317	-0.243234	0.292501	-0.263785	-0.194097	-0.337674	-0.119035	0.036482
4	4	0.643868	0.463769	0.912675	-0.136469	-0.955919	-0.282816	-0.506169	-0.169948	-0.242296	...	0.126639	-0.009587	0.236805	-0.441034	-0.283703	-0.013882	0.274577	-0.402516	-0.147622	0.223398

TPM034A Machine Learning for socio-technical systems¶

`Lab session 03: Working with image embeddings`¶

`Instructions`¶

`Use of AI tools`¶

`Google Colab workspace set-up`¶

`Application: Predicting perceived visual neighbourhood attractiveness`

`Organisation`¶

1. Loading and exploring the data¶

1.1. Reading the medatada and geospatial data¶

1.2. Exploring the image metadata¶

Exercise 1: Explore the statistics and geospatial distribution of the images¶

1.3. Exploring the geospatial data¶

Exercise 2: Do you recongnize some geospatial patterns of the attractiveness?¶

1.4. Visual inspection of images and attractiveness¶

2. Image embeddings¶

2.1. Embedding model¶

2.1. Exploring the embeddings¶

Exercise 3: Are year or month variables encoded by the images?¶

3. Training multiple models: Regression model, Random forest and Multilayer perceptron (MLP) for predicting attractiveness

3.1. Preparing the dataset¶

3.2. Linear multiple regression model¶

3.3. Random Forest (RF)¶

3.4 Multilayer perceptions (MLP)¶

Exercise 4: Apply GridSearch now using MLP to predict the attractiveness¶

3.5 Comparing and reflecting on the model performances and their outcomes¶

Exercise 5: Comparison of predictions¶

Exercise 6: Conclusion and reflection on the results¶

	img_id	Fhat000	Fhat001	Fhat002	Fhat003	Fhat004	Fhat005	Fhat006	Fhat007	Fhat008	...	Fhat015	Fhat016	Fhat017	Fhat018	Fhat019	Fhat020	Fhat021	Fhat022	Fhat023	Fhat024
count	101444.000000	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	...	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05	1.014440e+05
mean	50721.500000	-2.080693e-15	-9.357725e-16	-2.736434e-15	-4.974724e-15	2.683972e-15	7.304069e-16	1.346784e-15	-1.906707e-15	1.793097e-15	...	7.301267e-15	4.327248e-15	-1.516988e-15	-2.988028e-16	6.297553e-16	1.218466e-15	-2.973249e-15	9.047085e-16	-6.079300e-15	-3.081816e-15
std	29284.504691	6.683317e-01	6.462813e-01	5.985799e-01	5.460809e-01	5.299258e-01	4.759478e-01	4.310536e-01	3.771057e-01	3.599334e-01	...	2.749744e-01	2.690753e-01	2.675355e-01	2.605380e-01	2.513234e-01	2.370894e-01	2.330000e-01	2.261193e-01	2.154549e-01	2.138992e-01
min	0.000000	-2.205926e+00	-2.203099e+00	-2.196433e+00	-2.103170e+00	-1.533929e+00	-1.632734e+00	-1.761941e+00	-1.617488e+00	-1.396847e+00	...	-1.130731e+00	-1.117100e+00	-1.168965e+00	-1.157674e+00	-1.229036e+00	-1.202957e+00	-1.019485e+00	-1.166473e+00	-9.286276e-01	-1.018580e+00
25%	25360.750000	-4.672549e-01	-4.398960e-01	-4.065295e-01	-3.735691e-01	-3.803300e-01	-3.335537e-01	-2.915311e-01	-2.452553e-01	-2.427598e-01	...	-1.829394e-01	-1.765131e-01	-1.752905e-01	-1.691843e-01	-1.650595e-01	-1.546714e-01	-1.534472e-01	-1.424052e-01	-1.441687e-01	-1.389165e-01
50%	50721.500000	-3.127995e-02	-3.082165e-02	-1.821871e-03	2.867939e-03	-1.458150e-03	-3.181422e-02	-5.385117e-03	-6.570788e-03	-8.152784e-05	...	-4.343283e-03	-1.867820e-03	-8.741819e-03	-3.530588e-03	-3.853645e-03	1.124274e-03	-7.418046e-03	4.233806e-03	-6.926362e-03	-2.092290e-03
75%	76082.250000	4.443619e-01	4.109933e-01	4.118485e-01	3.631078e-01	3.636694e-01	3.016489e-01	2.858204e-01	2.359148e-01	2.437202e-01	...	1.770575e-01	1.747617e-01	1.656913e-01	1.645988e-01	1.618326e-01	1.560530e-01	1.436883e-01	1.458955e-01	1.359870e-01	1.372224e-01
max	101443.000000	2.526723e+00	3.405483e+00	2.218247e+00	2.476780e+00	2.073655e+00	2.156241e+00	1.930746e+00	1.950545e+00	1.422213e+00	...	1.617737e+00	1.422178e+00	1.513946e+00	1.329045e+00	1.243062e+00	1.530715e+00	1.342129e+00	1.246960e+00	1.083111e+00	1.164941e+00

TPM034A Machine Learning for socio-technical systems¶

Lab session 03: Working with image embeddings¶

Instructions¶

Use of AI tools¶

Google Colab workspace set-up¶

Application: Predicting perceived visual neighbourhood attractiveness

Organisation¶

1. Loading and exploring the data¶

1.1. Reading the medatada and geospatial data¶

1.2. Exploring the image metadata¶

Exercise 1: Explore the statistics and geospatial distribution of the images¶

1.3. Exploring the geospatial data¶

Exercise 2: Do you recongnize some geospatial patterns of the attractiveness?¶

1.4. Visual inspection of images and attractiveness¶

2. Image embeddings¶

2.1. Embedding model¶

2.1. Exploring the embeddings¶

Exercise 3: Are year or month variables encoded by the images?¶

3. Training multiple models: Regression model, Random forest and Multilayer perceptron (MLP) for predicting attractiveness

3.1. Preparing the dataset¶

3.2. Linear multiple regression model¶

3.3. Random Forest (RF)¶

3.4 Multilayer perceptions (MLP)¶

Exercise 4: Apply GridSearch now using MLP to predict the attractiveness¶

3.5 Comparing and reflecting on the model performances and their outcomes¶

Exercise 5: Comparison of predictions¶

Exercise 6: Conclusion and reflection on the results¶

`Lab session 03: Working with image embeddings`¶

`Instructions`¶

`Use of AI tools`¶

`Google Colab workspace set-up`¶

`Application: Predicting perceived visual neighbourhood attractiveness`

`Organisation`¶